Quasi-Static Fractures in Disordered Media and Iterated
Conformal Maps
Felipe Barra, H. George E. Hentschel*, Anders Levermann and Itamar Procaccia
Dept. of Chemical Physics, The Weizmann Institute of Science, Rehovot, 76100, Israel
* Dept of Physics, Emory University, Atlanta Georgia
Abstract
We study the geometrical cliaracteristic of quasi-static fractures in disordered media, using it-
erated conformal maps to determine the evolution of the fracture pattern. This method allows an
efficient and accurate solution of the Lame equations without resorting to lattice models. Typical
fracture patterns exhibit increased ramification due to the increase of the stress at the tips. We
find the roughness exponent of the experimentally relevant backbone of the fracture pattern; it
crosses over from about 0.5 for small scales to about 0.75 for large scales, in excellent agreement
with experiments. We propose that this cross-over reflects the increased ramification of the fracture
pattern.
Considerable amount of theoretical work |l], |], |^ on fracture in disordered media is based
on attempts to solve the equation of motion for an isotropic elastic body in the continuum
limit
= + '")^(^ • + -"^'^ • (1)
Here u is the field describing the displacement of each mass point from its location in an
unstrained body and p is the density. The constants fi and A are the Lame constants. In
terms of the displacement field the elastic strain tensor is defined as
'''-2 [dx, dxj ■
For the development of a crack the important object is the stress tensor, which in linear
elasticity is written as
aij = X6ij ^ ekk + 2fieij . (3)
k
When the stress component which is transverse to the interface of a crack exceeds a threshold
value (Tc, the crack can develop. When the external load is such that the transverse stress
exceeds only slightly the threshold value, the crack develops slowly, and one can neglect the
second time-derivative in Eq. (|l|). This is the quasi-static limit, in which after each growth
event one needs to recalculate the strain field by solving the Lame equation
(A + /i)V(V-u)+/iV2u = . (4)
In many previous works the problem was approached by discretizing Eq. on a lattice
@) i) S 01- ^^is letter we offer a novel approach based on iterated conformal maps; this
method turned out to be very useful in the context of fractal growth patterns p|, |^, [1^, |Tl
and it appears advantageous also for the present problem.
Although we can develop the approach in the full generality of Eq. (^, for the sake of
clarity in this Letter we will consider mode III fracturing for which a 3-dimensional elastic
medium is subjected to a finite shear stress (Xzy (y^o as y ±oo. Such an applied stress
will create a displacement field y), m^, = 0, Mj, = in the medium. Despite the medium
being three dimensional, therefore, the calculation of the strain and stress tensors are two
dimensional.
We can describe a crack of arbitrary shape by its interface x(s), where s is the arc length
which is used to parameterize the contour. We wish to develop a quasi-static model [12, 13
2
for the time development of this fracture in which discrete events advance the interface with
a normal velocity
Vn{s) = aia^tis) - crj (5)
if the transverse component of the stress tensor azt is greater than a critical yield value Uc for
fracturing; otherwise no fracture propagation occurs. We will use the notation {t, n) to de-
scribe respectively the transverse and normal directions at any point on the two-dimensional
crack interface. Whenever the interface has more than one position s for which Vn{s) does
not vanish, we will respect the disorder by choosing the next growth position randomly with
a probability proportional to Vn{s). This is similar to Diffusion Limited Aggregation (DLA)
in which a particle is grown with a probability proportional to the gradient to the field. One
should note that another model could be derived in which all eligible fracture sites are grown
simultaneously, growing a whole layer whose local width is Vn{s). This would be more akin
to Laplacian growth algorithms, which in general give rise to clusters in a different univer-
sality class than DLA [|^. We prefer the first in order to take the disorder of the material
explicitly into account.
In mode III fracture V ■ u = 0, and the Lame equation reduces to Laplace's equation
d^ujdx" + d^ujdy^ = 0, (6)
and therefore Uz is the real part of an analytic function
X{z) = Uz{x, y) + iiz{x, y) (7)
where z = x + iy. The boundary conditions far from the crack and on the crack interface
can be used to find this analytic function. It should be stated here that mode I and mode II
fractures can be reduced to a bi-Laplacian equation; then one needs to determine two rather
than one analytic functions. How to accomplish a growth model using iterated conformal
maps for those cases will be shown in a forthcoming publication |jl5[ .
Far from the crack as ?/ — ±oo we know cr^^ or using the stress/strain relationships
Eq. ^ we find that Uz ~ [aoo/ l^]y- Thus the analytic function must have the form
x{z) -i[a^/ii\z as \z\ ^ oo . (8)
Now on the boundary of the crack the normal stress vanishes, i.e.
= Oznis) = dnUz = -dt^z- (9)
Since is constant on the boundary, we choose ^2 = 0, which in turn is a boundary condition
making the analytic function xi^) real on the boundary of the crack:
X(^(.)) = xiz{s)r . (10)
The direct determination of the strain tensor for an arbitrary shaped (and evolving) crack
is still difficult. We therefore proceed by turning to a mathematical complex plane cu, in
which the crack is forever circular and of unit radius. The strain field for such a crack is
well known, being the real part of the function {u) where
X^'\u;)^-i[aJ^i](u;-l/u;) (11)
This is the unique analytic function obeying the boundary conditions x^°^(a;) —>■ —i[aoo/ fJ']uj
as — > 00, while on the unit circle x^^K^W^^) — X^^H^^P^^)*-
Now invoke a conformal map z — $(")(u;) that maps the exterior of unit circle in the
mathematical plane u to the exterior of the crack in the physical plane z, after n growth
steps. This conformal map is univalent by construction, and therefore admits a Laurent
expansion
$(")(a;) = Fi^'^Lj + F^^ + F^^^/uj + F^/^' + • • • (12)
Then the required analytic function x^'^\z) is given by the expression
X^"^^) = -i[F[^^<jJtJi\{^^-r\,) _ i/$W-\^)) (13)
Prom this we should compute now the transverse stress tensor:
dx^
n).
Z
ds
axW($W(e^^))a^
89 ds'
•pW d_( ie _ -i6\
|$'(n)(e^e)|
on the boundary.
Finally we describe how $^"■^(0;) is obtained. Suppose that ^'^'^~^\uj) is known, with
$*^°'(a;) being the identity, $(°^(a;) = u. We first compute the transverse strain tensor
4
o'ztid) = 2o"oo-^i"~"^'' cos^/|$'*^"~^-*(e~*^)|. In order to grow according to the requirement
we should choose growth sites more often when Acr(0) = cXztid) — (Jc is larger. We therefore
construct a probability density P{0) on the unit circle e*^ which satisfies
^ |$'(»-i)(e-^)|Aa(g)e(Aa(g))
where Q{Aa(6)) is the Heaviside function, and |$'("~^)(e*^)| is simply the Jacobian of the
transformation from mathematical to physical plane. The next growth position. On in the
mathematical plane, is chosen randomly with respect to the probability P{9)d9. In this way
the disorder of the medium is mirrored in the growth pattern. At the chosen position on
the crack, i.e. z = <|)("^i)(e*^"), we want to advance the crack with a region whose area
is the typical process zone for the material that we analyze. According to |Q the typical
scale of the process zone is K'^/a'^, where K is a characteristic fracture toughness parameter.
Denoting the typical area of the process zone by Aq, we achieve growth with an auxiliary
conformal map 0A„,e„(^^) that maps the unit circle to a unit circle with a bump of area A^
centered at e*^" . An example of such a map is given by |^ :
2w
X
l + w + w\ l + —
1 2 1-A\^/"
(Px^eiw) = e''<Pxfi{e~''w) , (17)
Here the bump has an aspect ratio a, < a < 1. In our work below we use a = 2/3. To
ensure a fixed size step in the physical domain we choose
A = ^0 (IQ)
|${"-i)'(e*^")|2 ■ ^
Finally the updated conformal map ^'^"^ is obtained as
=$(-l)(0,„^,J^)) . (19)
The recursive dynamics can be represented as iterations of the map 0a„,6»„(w^),
$(")(^) = <l>x^,e, o o . . . o 0A„,e„(^) . (20)
Every given fracture is determined completely by the random itinerary {6i}^^^. Eqs.(|T4D
together with (pO]) offer an analytic expression for the transverse stress field at any stage of
the crack propagation.
Fig. m exhibits a typical fracture pattern that is obtained with this theory, with cToo = 1,
after 10 000 growth events. The threshold value of cTc for the occurrence of the first event
(cf. Eq.(|T^) is (Tc = 2. We always implement the first event. For the next growth event
the threshold is (Tc = 2.9401.... We thus display in Fig.l a cluster obtained with (Jc = 2.94,
to be as close as possible to the quasi-static limit. Nevertheless, one should observe that
as the pattern develops, the stress at the active zone increases, and we get progressively
away from the quasi-static limit. Indeed, as a result of this, for fixed boundary conditions
at infinity, there are more and more values of 6 for which Eq. (|15D does not prohibit growth.
Since the tips of the patterns are mapped by ^ to larger and larger arcs on the unit
circle, the support of the probability P{6) increases, and the fracture pattern becomes more
and more ramified as the process advances. The geometric characteristics of the fracture
pattern are not invariant to the growth. For this reason it makes little sense to measure the
fractal dimension of the pattern; this is not a stable characteristic, and it will change with
the growth.
On the other hand, we should realize that the fracture pattern is not what is observed in
typical experiments. When the fracture hits the boundaries of the sample, and the sample
breaks into two parts, all the side-branches of the pattern remain hidden in the damaged
material, and only the backbone of the fracture pattern appears as the surface of the broken
parts. The backbone does not suffer from the geometric variability discussed above. In Fig.
^ we show the backbone of the pattern displayed in Fig. |I].
This backbone is representative of all the fracture patterns, we should note that in our
theory there are not lateral boundaries, and the backbone shown does not suffer from finite
size effects which may very well exist in experimental realizations.
In determining the roughness exponent of the backbone, we should note that a close
examination of it reveals that it is not a graph. There are overhangs in this backbone,
and since we deal with mode III fracturing, the two pieces of material can separate leaving
these overhangs intact. Accordingly, one should not approach the roughness exponent using
correlation function techniques; these may introduce serious errors when overhangs exist
T7|. Rather, we should measure, for any given r, the quantity |18
h{r) = (Max{?/(r')}x<r'<x+r - Min{y{r'}^,^r'<x+r)x ■ (21)
6
The roughness exponent ( is then obtained from
h{r) ~ , (22)
if this relation holds. To get good statistics we average, in addition to all x for the same
backbone, over many fracture patterns. The result of the analysis is shown in Fig. ^.
We find that the roughness exponent for the backbone exhibits a clear cross-over from
0.54 for shorter distances r to 0.75 for larger distances. Within the error bars these results
are in excellent agreement with the numbers quoted experimentally, see for example |]18
The short length scale exponent of order 0.5 is also in agreement with recent simulational
results of a lattice model [0 (which is by definition a short length scale solution). Bouchaud
||T8| proposed that the cross-over stems from transition between slow and rapid fracture, from
the "vicinity of the depinning transition" to the "moving phase" in her terms. Obviously,
in our theory we solve the quasi-static equation all along, and there is no change of physics.
Nevertheless, as we observed before, the fracture pattern begins with very low ramification
when the stress field exceeds the threshold value only at few positions on the fracture in-
terface. Later it evolves to a much more ramified pattern due to the increase of the stress
fields at the tips of the mature pattern. The scaling properties of the backbone reflect this
cross-over. We propose that this effect is responsible for the cross-over in the roughening
exponent of the backbone. It is remarkable that in spite of the apparently different mode of
fracture (mode III is difficult to realize experimentally) nevertheless the exponents appear
extremely close to those recorded experimentally for other fracture modes. This supports
the conjecture of universality proposed in
We have thus demonstrated that iterated conformal maps offer an efficient method for
studying fracture patterns. Here we considered only mode III quasi-static patterns. The
theory for mode I and mode II is available and will be presented elsewhere . The gener-
alization to dynamical scaling, in which Eq.(|l]) is considered including the time derivatives
is akin to the transition from electrostatics to electrodynamics. This is still an attractive
goal for the road ahead.
7
1000
y
-1000
-1000 1000
X
FIG. 1: A typical fracture pattern that is obtained from iterated conformal maps. What is seen is
the boundary of the fractured zone, which is the mapping of the unit circle in the mathematical
domain onto the physical domain. Notice that the pattern becomes more and more ramified as the
the fracture pattern develops. This is due to the enhancement of the stress field at the tips of the
growing pattern
Acknowledgments
We are indebted to S. Ciliberto for getting us interested in this problem, and to J.
Fineberg for some very useful discussions. This work has been supported in part by the
Petroleum Research Fund, The European Commission under the TMR program and the
Naftali and Anna Backenroth-Bronicki Fund for Research in Chaos and Complexity. A. L.
is supported by a fellowship of the Minerva Foundation, Munich, Germany.
[1] N.I. Muskhelishvili, Some Basic Problems in the Mathematical Theory of Elasticity, (Noord-
hoff, Groningen, 1952).
[2] L.D. Landau and E.M. Lifshitz, Theory of Elasticity, 3rd ed. (Pergamon, London, 1986).
[3] J. Fineberg and M. Marder, Phys. Rep. 313, 1, (1999), and references therein.
[4] H.J. Herrmann and S. Roux, Statistical Models for the Fracture of Disordered Media, (North
Holland, Amsterdam, 1990), and references therein.
[5] L. Slepian, Dov. Phys. Dokl. 26, 538 (1981).
8
1000
y
-1000
-1000
1000
FIG. 2: A typical backbone of the fracture pattern. This is the projection onto the x-y plane of
the experimentally observed boundary between the two parts of the material that separate when
the fracture pattern hits t
1000
FIG. 3: /i(r) averaged over all the backbone and over 70 fracture patterns each of which of 10 000
fracture events. There is a cross-over between a scaling law with roughness exponent 0.54 it 0.05
to and exponent of 0.75 ± 0.02
[6] M. Marder and X. Liu, Phys. Rev. Lett. 71, 2417 (1993).
[7] A. Parisi, G. Caldarelli and L. Pietronero, arXiv:cond-mat / 0004374
[8] M.B. Hastings and L.S. Levitov, Physica D 116, 244 (1998).
[9] B. Davidovitch, H.G.E. Hentschel, Z. Olami, I.Procaccia, L.M. Sander, and E. Somfai, Phys.
Rev. E, 59 1368 (1999).
[10] B. Davidovitch, M.J. Feigenbaum, H.G.E. Hentschel and I. Procaccia, Phys. Rev. E 62, 1706
(2000).
[11] B. Davidovitch, A. Levermann, I. Procaccia, Phys. Rev. E, 62 R5919.
[12] M. Barber, J. Donley and J. S. Langer, Phys. Rev. A40, 366 (1989)
[13] See for example J. Kertesz in |Q.
[14] F. Barra, B. Davidovitch and I. Procaccia, "Iterated Conformal Dynamics and Laplacian
Growth", submitted to Phys. Rev. E, also cond-mat /0105608
[15] F. Barra, A. Levermann and I. Procaccia, "Solution of the bi-Laplacian Equation for Quasi-
Static Mode I and Mode II Fracture Using Iterated Conformal Maps" , in preparation.
[16] See for example [^ Sect. 2.3.
[17] Z. Olami, I. Procaccia and R. Zeitak, Phys. Rev. E52, 3402 (1995).
[18] E. Bouchaud, J. Phys: Condens. Matter 9, 4319 (1997).
10